Python para Otimização - SCIP

Leitura e escrita em arquivos

A primeira coisa que precisamos fazer é conseguir ler dados de um arquivo externo, pois serão a partir desses dados que passaremos os inputs dos nossos problemas de otimização.

Introdução

Vamos aprender um método que lê os dados de um arquivo linha-a-linha. Considere um arquivo de texto com 4 linhas, da seguinte forma:

Podemos o usar o seguinte código para ler a primeira linha:

# Caminho do arquivo (não esquecer o 'r' antes de começar)
caminho_arquivo = r"G:\Meu Drive\Arquivos\UFPR\1 - Disciplinas\8 - PI Avançada\Codigos\Python\1-LeituraDeDados\input.txt"

# Abre o arquivo
with open(caminho_arquivo, "r", encoding="utf-8") as arquivo:
    # Atribui a primeira linha na variável 'a'
    a = arquivo.readline()
    # imprime a variável
    print(a)
## Linha 1

A cada vez que executamos o comando arquivo.readline() uma nova linha é lida. Como o arquivo possui 4 linhas, no exemplo abaixo executamos e imprimimos 4 vezes:

caminho_arquivo = r"G:\Meu Drive\Arquivos\UFPR\1 - Disciplinas\8 - PI Avançada\Codigos\Python\1-LeituraDeDados\input.txt"
with open(caminho_arquivo, "r", encoding="utf-8") as arquivo:
    a = arquivo.readline()
    print(a)
    a = arquivo.readline()
    print(a)
    a = arquivo.readline()
    print(a)
    a = arquivo.readline()
    print(a)
## Linha 1
## 
## 10
## 
## Linha 3
## 
## 2 4

Ou ainda, podemos imprimir usando um laço for

caminho_arquivo = r"G:\Meu Drive\Arquivos\UFPR\1 - Disciplinas\8 - PI Avançada\Codigos\Python\1-LeituraDeDados\input.txt"
with open(caminho_arquivo, "r", encoding="utf-8") as arquivo:
    for i in range(4):
      a = arquivo.readline()
      print(a)
## Linha 1
## 
## 10
## 
## Linha 3
## 
## 2 4

Usando indicadores no arquivo

Muitas vezes, para criar programas que são genéricos, não sabemos com antecedência o número de linhas de um arquivo, e portanto, precisamos indicar isso no próprio arquivo. Considere a leitura de vários números que estão em um arquivo, porém, o número de elementos a serem lidos é indicado no próprio arquivo, como seu primeiro elemento. Nesse caso, precisamos criar uma lógica que lê o primeiro elemento e depois cria o laço para ler o restante. Considere um arquivo em que precisamos ler os números 10, 20, 30, 40 e 50. Podemos usar o seguinte arquivo para tal (note que são 5 números, portanto o primeiro elemento do arquivo é 5):

caminho_arquivo = r"G:\Meu Drive\Arquivos\UFPR\1 - Disciplinas\8 - PI Avançada\Codigos\Python\1-LeituraDeDados\input2.txt"
with open(caminho_arquivo, "r", encoding="utf-8") as arquivo:
    n = arquivo.readline() # lê o primeiro elemento (como string)
    n = int(n)             # transforma o elemento em númerico (inteiro)
    for i in range(n):
      a = arquivo.readline()
      print(a)
## 10
## 
## 20
## 
## 30
## 
## 40
## 
## 50

Observe que foi necessário transformar o dado em inteiro para usar na função range(). Isso ocorre pois sempre que lêmos um dado de um arquivo externo, ele é armazenado como uma string.

Lendo múltiplos elementos em uma linha

Até agora fizemos a leitura linha-a-linha, sendo que cada elemento ocupava ma linha. Mas e se existirem vários elementos em uma mesma linha? Considere o seguinte arquivo, em que devemos ler as palavras “Esta é uma frase”, e todas estão na mesma linha:

caminho_arquivo = r"G:\Meu Drive\Arquivos\UFPR\1 - Disciplinas\8 - PI Avançada\Codigos\Python\1-LeituraDeDados\input3.txt"
with open(caminho_arquivo, "r", encoding="utf-8") as arquivo:
    frase = arquivo.readline() 
    print(frase)
## Está é uma frase

Mas que tipo de estrutra de dados tem a variável frase, e como poderíamos acessar os elementos de forma individual? A função type mostra o tipo de dados de uma estrutura:

caminho_arquivo = r"G:\Meu Drive\Arquivos\UFPR\1 - Disciplinas\8 - PI Avançada\Codigos\Python\1-LeituraDeDados\input3.txt"
with open(caminho_arquivo, "r", encoding="utf-8") as arquivo:
    frase = arquivo.readline() 
    print(type(frase))
## <class 'str'>

Como anteriormente, toda a frase foi lida em uma string. Podemos acessar os elementos indivíduais usando a função split() da própria string. A função retorna uma list, que cada elemento é uma substring, e é possível definir qual o separador usado (espaço, ponto-e-vírgula, etc…).

O código abaixo mostra a separação da primeira linha na lista l_frase, e em seguida todos os elementos são impressos usando um laço for:

caminho_arquivo = r"G:\Meu Drive\Arquivos\UFPR\1 - Disciplinas\8 - PI Avançada\Codigos\Python\1-LeituraDeDados\input3.txt"
with open(caminho_arquivo, "r", encoding="utf-8") as arquivo:
    frase   = arquivo.readline() 
    l_frase = frase.split(" ") # Usando split com separador " " (espaço)
    # Imprimindo os elementos, um a um
    for e in l_frase:
      print(e)
## Está
## é
## uma
## frase

OBS: Para saber mais sobre listas, veja o material, na seção 3.2.2 - Listas.

Escrevendo em arquivos

Além de ler dados de arquivos de texto, também precisamos entender como escrever. Quando resolvemos um modelo de otimização, podemos usar os arquivos de texto para exportar as soluções. O código abaixo abre um novo arquivo e escreve uma frase nele:

caminho_arquivo = r"G:\Meu Drive\Arquivos\UFPR\1 - Disciplinas\8 - PI Avançada\Codigos\Python\1-LeituraDeDados\output.txt"
with open(caminho_arquivo, "w", encoding="utf-8") as arquivo:
  arquivo.write("Frase escrita no arquivo, direto do python")

Para pular linhas escrevemos a string \n:

caminho_arquivo = r"G:\Meu Drive\Arquivos\UFPR\1 - Disciplinas\8 - PI Avançada\Codigos\Python\1-LeituraDeDados\output2.txt"
with open(caminho_arquivo, "w", encoding="utf-8") as arquivo:
  arquivo.write("Frase\n escrita\n em\n várias\n linhas\n")

Instalação e Introdução SCIP

Para instalar o pacote do scip (pyscipopt usando o pip, basta usar o comando no terminal:

pip install pyscipopt

Criando um modelo

O principal objeto do solver é o modelo, que é feito da seguinte forma:

# Importando o modelo
from pyscipopt import Model

# Criando um objeto do tipo Model()
scip = Model()

Agora podemos popular um modelo, adicionando varíaveis e restrições. Vamos implementar o modelo abaixo (problema do sapateiro):

\[\begin{alignat*}{4} & \text{max z = } & 5x_1 & + 2x_2 \\ & \text{Sujeito à} & 10x_1 & + 12x_2 &\leq 60\\ & & 2x_1 & + x_2 & \leq 6\\ & & x_1 \geq 0 & \quad x_2 \geq 0& \end{alignat*}\]

from pyscipopt import Model

# Criando um objeto do tipo Model()
scip = Model()

# Declarando as variáveis
x1 = scip.addVar(vtype = 'C', lb = 0, ub = None, name = 'x1')
x2 = scip.addVar(vtype = 'C', lb = 0, ub = None, name = 'x2')

# Restrições
res1 = scip.addCons(10*x1 + 12*x2 <= 60, name = "res1")
res2 = scip.addCons(2*x1  + x2    <= 6,  name = "res2")

# Função objetivo
scip.setObjective(5*x1 + 2*x2, sense = "maximize")

# Otimizando o modelo
scip.optimize()

vtype = "C" indica uma variável contínua, poderia ser usado "I" (inteira) e "B" (binária). lb e ub indicam os limites inferiores e superiores das variáveis, neste caso de 0 a infinito.

Coletando os resultados

Agora podemos coletar os valores de função objetivo e variáveis:

z      = scip.getObjVal()
x1_val = scip.getVal(x1)
x2_val = scip.getVal(x2)
print("Valor objetivo :{}\nx1: {}\nx2: {}".format(z,x1,x2))
## Valor objetivo :15.0
## x1: x1
## x2: x2

Exercícios

Variáveis

No exemplo acima, criamos duas variáveis do Pyhton para armazenar as duas variáveis do SCIP (x1 e x2). Para modelos maiores, no entanto, não é viável definir as variáveis dessa forma, sendo que usamos as próprias estruturas de dados do Python para definiar conjuntos maiores de variáveis dos modelos.

Vamos considerar o seguinte problema do transporte:

(O problema do transporte) Uma empresa possui 2 fábricas (I e II) e 3 depósitos (A,B e C). Cada fábrica possui uma capacidade de produção e cada depósito uma demanda de consumo. Existe uma distância entre cada fábrica e cada depósito, de forma que há um custo associado a cada unidade de produto transportado de uma fábrica a um depósito. A empresa precisa atender às demandas dos depósitos, sem exceder as capacidades produtivas das fábricas. Os dados de demandas, capacidades e custos de transporte são mostrados na Tabela abaixo:

Vamos considerar que teremos variáveis \(x_{00}, x_{01}, x_{02}, x_{10}, x_{11}, x_{12}\).

Dicionários

from pyscipopt import Model
scip = Model()

var_dict    = {} 
n_fabricas  = 2
n_depositos = 3

# Como existem dois niveis de variaveis, criamos um dicionário de dicionários
for i in range(n_fabricas):
  var_dict[i] = {}
  for j in range(n_depositos):
    var_dict[i][j] =  scip.addVar(vtype = 'C', name = f"x_{i}_{j}")

print(var_dict)
## {0: {0: x_0_0, 1: x_0_1, 2: x_0_2}, 1: {0: x_1_0, 1: x_1_1, 2: x_1_2}}

scip.freeProb()
scip = None

Listas

Da mesma forma, podemos criar uma lista de listas para armazenar as variáveis:

from pyscipopt import Model
scip = Model()

n_fabricas  = 2
n_depositos = 3
var_list    = [ [None for j in range(n_depositos)] for i in range(n_fabricas)] 

for i in range(len(var_list)):
  for j in range(len(var_list[0])):
    var_list[i][j] =  scip.addVar(vtype = 'C', name = f"x_{i}_{j}")

print(var_list)
## [[x_0_0, x_0_1, x_0_2], [x_1_0, x_1_1, x_1_2]]

scip.freeProb()

Restrições e Quicksum

Da mesma forma que podemos criar uma lista com centenas de variáveis, podemos criar restrições com essas variáveis no formato de somatório usando o quicksum. Considere que as restrições de capacidade das fábricas devem ficar da seguinte forma:

\[\begin{alignat*}{4} & & x_{00} & + x_{01} & + x_{02} \leq 350\\ & & x_{10} & + x_{11} & + x_{12} \leq 650\\ \end{alignat*}\]

from pyscipopt import Model, quicksum
scip = Model()

# Restrição de demanda 1
scip.addCons(  quicksum(var_list[0][j] for j in range(n_depositos)) <= 350 , name = "res1")
## res1
# Restrição de demanda 1
scip.addCons(  quicksum(var_list[1][j] for j in range(n_depositos)) <= 650 , name = "res2")
## res2

Também podemos usar o quicksum para criar a função objetivo. Abaixo o modelo do transporte completo:

from pyscipopt import Model, quicksum
scip = Model()


ll_custos      = [[2.5,1.7,1.8],[2.5,1.8,1.4]]
l_demandas     = [300,300,300]
l_capacidades  = [350,650]
n_fabricas  = 2
n_depositos = 3

# Variáveis
var_list    = [ [None for j in range(n_depositos)] for i in range(n_fabricas)] 
for i in range(len(var_list)):
  for j in range(len(var_list[0])):
    var_list[i][j] =  scip.addVar(vtype = 'C', name = f"x_{i}_{j}")

# Restrições de capacidade
for i in range(n_fabricas):
  scip.addCons(quicksum(var_list[i][j] for j in range(n_depositos)) <= l_capacidades[i], name = "cap. {}".format(i))
  
# Restrições de demanda
for j in range(n_depositos):
  scip.addCons(quicksum(var_list[i][j] for i in range(n_fabricas))  >= l_demandas[i], name = "dem. {}".format(j))

# Função objetivo
scip.setObjective( quicksum(ll_custos[i][j] * var_list[i][j] for i in range(n_fabricas) for j in range(n_depositos)), sense = "minimize")

scip.optimize()
z     = scip.getObjVal()
print(z)

Exportando modelos

No scip, podemos exportar modelos em diversos formatos (.lp, .mps, .sol, etc… https://pyscipopt.readthedocs.io/en/stable/tutorials/readwrite.html). Para fazer a exportação usamos a função writeProblem:

from pyscipopt import Model
scip = Model()
scip.writeProblem(filename = "example_file.mps", trans = False, genericnames = False)
#mestre.writeProblem(r"G:\\Meu Drive\\Arquivos\\UFPR\\1 - Disciplinas\\8 - PI Avançada\\Codigos\\Python\\1-LeituraDeDados\\modelo.lp")

Eu já tive problemas com alguns caminhos na exportação, consegui resolver quando coloquei a letra r antes da string, e usei barras duplas, como no comentário do código acima.

Dualidade

Muitos métodos avançados precisam extrair os valores duais, e a maioria dos solvers fornece essa funcionalidade. No SCIP, no entanto, existem alguns “problemas” inerentes à coleta dos valores duais, em função da própria forma como ele foi implementado (alguns lugares que tratam disso: aqui e aqui).

Básicamente, para conseguirmos extrair os valores duais do SCIP um conjunto de configurações deve ser desabilitado (aqueles que permitem ao SCIP realizar tranformações no modelo, para fins de eficiência computacional). São eles:

from pyscipopt import Model, SCIP_PARAMSETTING
scip = Model()
scip.setPresolve(SCIP_PARAMSETTING.OFF)
scip.setHeuristics(SCIP_PARAMSETTING.OFF)
scip.disablePropagation()

Associado a cada restrição de um PL temos um valor dual, que pode ser coletado usando a função getDualsolLinear(rest), em que rest é a restrição. Considere novamente o modelo do sapateiro, apresentado na seção “Criando um modelo”.

from pyscipopt import Model, SCIP_PARAMSETTING

# Criando um objeto do tipo Model()
scip = Model()

scip.setPresolve(SCIP_PARAMSETTING.OFF)
scip.setHeuristics(SCIP_PARAMSETTING.OFF)
scip.disablePropagation()

# Declarando as variaveis
x1 = scip.addVar(vtype = 'C', lb = 0, ub = None, name = 'x1')
x2 = scip.addVar(vtype = 'C', lb = 0, ub = None, name = 'x2')

# Restricoes
res1 = scip.addCons(10*x1 + 12*x2 <= 60, name = "res1")
res2 = scip.addCons(2*x1  + x2    <= 6,  name = "res2")

# Funcao objetivo
scip.setObjective(5*x1 + 2*x2, sense = "maximize")

# Otimizando o modelo
scip.optimize()

# Coletando os valores duais:
pi1 = scip.getDualsolLinear(res1)
pi2 = scip.getDualsolLinear(res2)

print("pis",pi1,pi2)

CUIDADO: Existe ainda uma situação em que os valores duais podem estar errados - em restrições com somente uma variável. Nesses casos, ao que parece, mesmo informando que o SCIP não deve fazer transformações no modelo, ele o faz mesmo assim, considerando a restrição como uma variável limitada, e simplificando a restrição. Uma saída para isso é simplesmente criar uma variável “Dummy”, com lb = ub = 0, e adicionar ela nas restrições.

Considere os valores do exemplo abaixo com o problema do sapateiro alterando a segunda restrição, para que ela possua somente uma variável:

from pyscipopt import Model, SCIP_PARAMSETTING

# Criando um objeto do tipo Model()
scip = Model()

scip.setPresolve(SCIP_PARAMSETTING.OFF)
scip.setHeuristics(SCIP_PARAMSETTING.OFF)
scip.disablePropagation()

# Declarando as variaveis - com dummy
x1     = scip.addVar(vtype = 'C', lb = 0, ub = None, name = 'x1')
x2     = scip.addVar(vtype = 'C', lb = 0, ub = None, name = 'x2')
dummy  = scip.addVar(vtype = 'C', lb = 0, ub = 0, name = 'dummy')
# Restricoes
res1 = scip.addCons(10*x1 + 12*x2        <= 60, name = "res1")
res2 = scip.addCons(2*x1  + dummy        <= 6,  name = "res2")

# Funcao objetivo
scip.setObjective(5*x1 + 2*x2, sense = "maximize")

# Otimizando o modelo
scip.optimize()

# Coletando os valores duais:
pi1 = scip.getDualsolLinear(res1)
pi2 = scip.getDualsolLinear(res2)

print("pis",pi1,pi2)

E agora para corrigir, adicionamos uma variável dummy nessa restrição, e os valores duais ficam corretos.

Tópicos avançados

Decomposição de Dantzig-Wolfe

Exemplo Simples

import logging

logger = logging.getLogger(__name__)
logger.setLevel(logging.DEBUG)

console_handler = logging.StreamHandler()
console_handler.setLevel(logging.DEBUG)

class LoggerColor(logging.Formatter):
    def format(self, record):
        msg = super().format(record)

        if record.levelname == "DEBUG":
            levelname = f"\033[94m\033[1m==> {record.levelname}:\033[0m" # blue
        elif record.levelname == "INFO":
            levelname = f"\033[92m\033[1m=> {record.levelname}:\033[0m" # green
        elif record.levelname == "WARNING":
            levelname = f"\033[93m\033[1m===> {record.levelname}:\033[0m"  # yellow
        elif record.levelname == "ERROR":
            levelname = f"\033[38;5;214m\033[1m===> {record.levelname}:\033[0m"  # orange
        elif record.levelname == "CRITICAL":
            levelname = f"\033[91m\033[1m===> {record.levelname}:\033[0m" # red
        else:
            levelname = record.levelname 

        return f"{levelname} {record.message}"

formatter = LoggerColor("%(levelname)s: %(message)s")
console_handler.setFormatter(formatter)

logger.addHandler(console_handler)
# Exemplo de aplicação da decomposição de Dantzig-Wolfe
# usando o SCIP

import logging

from numpy import pi
from logger_config import logger
from pyscipopt import Model, quicksum, SCIP_PARAMSETTING

def gera_col(sol_x, n_cons):
    c    = 0
    x    = [None for i in range(n_cons)]
    c    = sol_x[0] - 2 * sol_x[1]
    x[0] =   2 * sol_x[0] + 2 * sol_x[1]
    x[1] = - 2 * sol_x[0] + 2 * sol_x[1]
    x[2] =   2 * sol_x[0] + 1 * sol_x[1]
    x[3] =   1

    return c, x

def calcula_fo_subproblema(pi, n_var):
    custo = [None for i in range(n_var)]
    custo[0] =  1 - 2 * pi[0] + 2 * pi[1] - 2 * pi[2] # x1
    custo[1] = -2 - 2 * pi[0] - 2 * pi[1] - 1 * pi[2] # x2
    return custo

def gera_mestre_initial_col(inital_sol, initial_c):
    mestre = Model("Mestre Dantzig-Wolfe")
    mestre.hideOutput(True)
    #mestre.enableReoptimization()
    mestre.setPresolve(SCIP_PARAMSETTING.OFF)
    mestre.setHeuristics(SCIP_PARAMSETTING.OFF)
    mestre.disablePropagation()

    # A última é a var. Dummy
    mestre_dummy_var = mestre.addVar("dummy", vtype = "C", lb = 0)
    # Inicializa com uma coluna (uma variável)
    alpha            = [mestre.addVar("a0", vtype = "C", lb = 0)]
    sol              = inital_sol # Initial sol => x1 x2
    c, nova_coluna   = gera_col(sol, 4)
    #logger.info(f"Custo da nova coluna:, {c}")
    #logger.info(f"Nova coluna:, {nova_coluna}")

    # Adiciona a primeira coluna ao mestre
    # Custo da fo
    mestre.setObjective(c * alpha[0], sense = "minimize")

    # Restrições do mestre
    mestre.addCons(alpha[0] * nova_coluna[0] + mestre_dummy_var >= 3)
    mestre.addCons(alpha[0] * nova_coluna[1] + mestre_dummy_var <= 3)
    mestre.addCons(alpha[0] * nova_coluna[2] + mestre_dummy_var <= 10)
    mestre.addCons(alpha[0] * nova_coluna[3] + mestre_dummy_var == 1)
    return mestre, alpha, mestre_dummy_var

def gera_subproblema_inicial():
    subproblema = Model("Sub Dantzig-Wolfe")
    subproblema.hideOutput(True)
    subproblema.setPresolve(SCIP_PARAMSETTING.OFF)
    subproblema.setHeuristics(SCIP_PARAMSETTING.OFF)
    subproblema.disablePropagation()

    # Variáveis do subproblema
    x = [subproblema.addVar(f"x_{i}" ,vtype = "C", lb = 0) for i in range(2)]
    # Restrições do subproblema
    subproblema.addCons(1 <= (x[0] <= 3))
    subproblema.addCons(1 <= (x[1] <= 3))
    # Fo
    subproblema.setObjective(0, sense = "minimize")
    return subproblema, x

def add_col_mestre(mestre, alpha, c, nova_coluna):
    # Add. a variável e o custo na fo
    alpha.append(mestre.addVar(f"a{len(alpha)}", vtype = "C", lb = 0, obj = c))
    i = 0
    for c in mestre.getConss():
        mestre.addConsCoeff(c, alpha[-1], nova_coluna[i])  # adds 1.0 * z to cons1
        i = i + 1
    return mestre,alpha

# Mestre inicializado com uma sol. inicial
initial_sol = [2,2]
mestre, alpha, mestre_dummy_var = gera_mestre_initial_col(initial_sol, 0)

# Subproblema inicial com fo sem termos
sub, x = gera_subproblema_inicial()
l_pontos = []
l_pontos.append(initial_sol)
custo_atualizado = -float('inf')
i = 0
while custo_atualizado + 0.01 < 0:
    mestre.writeProblem(f"C:\\Users\\x-eco\\OneDrive - ufpr.br\\Documentos\\Alexandre\\DW\\mestre_{i}.lp",trans = False, genericnames = False)
    mestre.optimize()
    # Duais, fo e variáveis
    pi          = [mestre.getDualsolLinear(c) for c in mestre.getConss()] # keep dual variables
    z           = mestre.getObjVal()
    alpha_val   = [mestre.getVal(v) for v in alpha]
    logger.info(f"Sol. Ótima do Mestre: {alpha_val}")
    logger.info(f"Custo: {z}")
    logger.info(f"Valores duais: {pi}")

    # Calcula novo custo do subproblema
    custo_sub = calcula_fo_subproblema(pi, len(x))
    sub.freeTransform()
    sub.chgReoptObjective(quicksum(custo_sub[i] * x[i] for i in range(len(x))), sense = 'minimize')
    sub.writeProblem(f"C:\\Users\\x-eco\\OneDrive - ufpr.br\\Documentos\\Alexandre\\DW\\sub_{i}.lp",trans = False, genericnames = False)
    sub.optimize()

    x_sol = [sub.getVal(x[i]) for i in range(len(x))]
    z_sub = sub.getObjVal()
    logger.info(f"Sol. Ótima do Subproblema: {x_sol}, com valor: {z_sub}")
    l_pontos.append(x_sol)

    custo_atualizado = (z_sub - pi[-1])
    logger.warning(f"Custo atualizado: {custo_atualizado}")
    if custo_atualizado > -0.01:
        logger.info("Solução ótima encontrada pelo Mestre")
    else:
        mestre.freeTransform()
        logger.info("Coluna consegue melhorar a solução do Mestre, adicionando nova coluna...")
        c, nova_coluna = gera_col(x_sol, 4)
        mestre,alpha = add_col_mestre(mestre, alpha, c, nova_coluna)
        #mestre.writeProblem(r"G:\\Meu Drive\\Arquivos\\UFPR\\1 - Disciplinas\\8 - PI Avançada\\Codigos\\Python\\1-LeituraDeDados\\modelo.lp")
    i = i+1
logger.info("Processo de Dantzig-Wolfe finalizado!!")
logger.info(f"Pontos gerados: {l_pontos}")

Exemplo com estrutura em blocos

# Exemplo de decomposição de DW usando
# a estrutura em blocos


# Exemplo de aplicação da decomposição de Dantzig-Wolfe
# usando o SCIP

import logging

from logger_config import logger
from pyscipopt import Model, quicksum, SCIP_PARAMSETTING


def gera_mestre_initial_col(param, initial_sol_x, initial_sol_I):
    
    mestre = Model("Mestre Dantzig-Wolfe")
    #mestre.hideOutput(True)
    mestre.setPresolve(SCIP_PARAMSETTING.OFF)
    mestre.setHeuristics(SCIP_PARAMSETTING.OFF)
    mestre.disablePropagation()

    """"
    alpha            = [ [mestre.addVar(f"a_{i}_0", vtype = "C", lb = 0)] for i in range(len(param["D"])) ]
    var_art          = mestre.addVar(f"va", vtype = "C", lb = 0)
    mestre_dummy_var = mestre.addVar("dummy", vtype = "C", lb = 0, ub = 0)
    
    mestre.addCons(var_art + 450 * alpha[0][0] + 32 * alpha[1][0] + mestre_dummy_var <= 240)
    mestre.addCons(var_art +  112 * alpha[1][0] + mestre_dummy_var <= 320) 
    mestre.addCons(var_art  + mestre_dummy_var <= 200)
    mestre.addCons(var_art + alpha[0][0] + mestre_dummy_var == 1)
    mestre.addCons(var_art + alpha[1][0] + mestre_dummy_var == 1)

    mestre.setObjective(100000 * var_art + 6750 * alpha[0][0] + 1100 * alpha[1][0], sense = "minimize")
    """
    
    n_prod = len(param["D"]) # número de produtos
    var_art          = mestre.addVar(f"va", vtype = "C", lb = 0)
    # Inicializa uma variável para cada produto
    alpha           = [ [mestre.addVar(f"a_{i}_0", vtype = "C", lb = 0)] for i in range(len(param["D"])) ]
    
    # A última é a var. Dummy
    mestre_dummy_var = mestre.addVar("dummy", vtype = "C", lb = 0, ub = 0)

    # Adicionando uma coluna por produto
    l_c      = []
    l_coluna = []
    for i in range(n_prod):
        # Cria a coluna por produto i
        c, nova_coluna  = gera_col(param, initial_sol_x[i], initial_sol_I[i], i)
        l_c.append(c)
        l_coluna.append(nova_coluna)

    # Adiciona restrições somente com a var alpha 1 (para ter um RHS), depois adiciona por inserção de coluna mesmo
    # primeiro as restrições de capacidade dos tanques
    for j in range(len(param["DispTanque"])):
        mestre.addCons(alpha[0][0] * l_coluna[0][j] + mestre_dummy_var <= param["DispTanque"][j])

    # Agora as restrições de convexidade
    for j in range(len(param["DispTanque"]),len(param["DispTanque"]) + len(param["D"])): # Convexidade
        mestre.addCons(alpha[0][0] * l_coluna[0][j] + mestre_dummy_var == 1)

    # Adiciona as outras colunas por inserção de coluna
    for p in range(1, len(param["D"])): # para cada produto além do primeiro
        for i,c in enumerate(mestre.getConss()):     # para cada restrição
            mestre.addConsCoeff(c, alpha[p][0], l_coluna[p][i])
 
    # Custo da fo
    mestre.setObjective(quicksum(l_c[i] * alpha[i][0] for i in range(n_prod)), sense = "minimize")
    
    return mestre, alpha, var_art

def gera_col(param, sol_x, sol_I, sub_prob):
    
    """
    Gera uma nova coluna a partir da solução do subproblema sol_x e sol_I
    sub_prob indica a qual subproblema a coluna pertence, pois isso vai influenciar
    na posição em que o valor da restrição de convexidade é adicionado.
    Finalmente, n_cons indica o número total de restrições do mestre.
    """

    n_prod = len(param["D"]) # número de produtos

    c_coeff       = 0
    alpha_coeff   = [0 for i in range(len(param["DispTanque"]) + n_prod) ]

    for i in range(len(param["DispTanque"])):
        alpha_coeff[i] = param["HMaquina"][sub_prob] * sol_x[i]

    alpha_coeff[len(param["DispTanque"]) + sub_prob] = 1  # Restrição de convexidade

    c_coeff       = sum(x * c for x, c in zip(sol_x, param["Cp"][sub_prob])) 
    c_coeff       = c_coeff + sum(I * c for I, c in zip(sol_I, param["Ci"][sub_prob]))
    return c_coeff, alpha_coeff

def gera_subproblem(param, sub_prob):

    n_prod = len(param["D"]) # número de produtos
    T = len(param["D"][0])   # número de períodos

    subproblema = Model(f"sub_{sub_prob}")
    subproblema.hideOutput(True)
    subproblema.setPresolve(SCIP_PARAMSETTING.OFF)
    subproblema.setHeuristics(SCIP_PARAMSETTING.OFF)
    subproblema.disablePropagation()

    # Variáveis do subproblema
    x = [subproblema.addVar(f"x_{sub_prob}_{i}" ,vtype = "C", lb = 0) for i in range(T)]
    I = [subproblema.addVar(f"I_{sub_prob}_{i}" ,vtype = "C", lb = 0) for i in range(T)]
    
    # Restrições do subproblema
    for t in range(T):
        demanda_t = param["D"][sub_prob][t]
        if t == 0:
            subproblema.addCons(x[t] - I[t] == demanda_t)
        else:
            subproblema.addCons(x[t] + I[t-1] - I[t] == demanda_t)

    # Fo
    subproblema.setObjective(0, sense = "minimize")
    return subproblema, x, I

def add_col_mestre(mestre, alpha, c, nova_coluna, sub_prob):
    # Add. a variável e o custo na fo
    alpha[sub_prob].append(mestre.addVar(f"a_{sub_prob}_{len(alpha[sub_prob])}", vtype = "C", lb = 0, obj = c))
    i = 0
    for i,c in enumerate(mestre.getConss()):
        mestre.addConsCoeff(c, alpha[sub_prob][-1], nova_coluna[i])  # adds 1.0 * z to cons1
    return mestre, alpha

def recupera_solucao(alpha_val, l_solucoes_x, l_solucoes_I):
    n_prod = len(l_solucoes_x)
    T      = len(l_solucoes_x[0][0])
    solucoes_x = [ [0 for t in range(T)] for p in range(n_prod)]
    solucoes_I = [ [0 for t in range(T)] for p in range(n_prod)]

    for p in range(n_prod):
        for t in range(T):
            for k in range(len(l_solucoes_x[p])): # para cada coluna usada
                solucoes_x[p][t] += alpha_val[p][0][k] * l_solucoes_x[p][k][t]
                solucoes_I[p][t] += alpha_val[p][0][k] * l_solucoes_I[p][k][t]
    return solucoes_x, solucoes_I

param = {"D":[[900,1800,1800],[400,600,800]],
        "Cp":[[1.0,1.5,2.0],[0.5,0.5,0.9]],
        "Ci":[[0.5,0.25,0],[0.25,0.25,0]],
        "DispTanque":[240,320,200],
        "HMaquina":[0.1,0.08],}

n_prod = len(param["D"]) # número de produtos
T      = len(param["D"][0])   # número de períodos

# Mestre inicializado com uma sol. inicial
initial_sol_x = [[1800,1800,900],[750,250,800]]
initial_sol_I = [[900,900,0],[350,0,0]]

mestre, alpha, var_art = gera_mestre_initial_col(param, initial_sol_x, initial_sol_I)
mestre.writeProblem(f"C:\\Users\\x-eco\\OneDrive - ufpr.br\\Documentos\\Alexandre\\DW\\mestre_0.lp", trans = False, genericnames = False)
alpha_val   = [[] for i in range(n_prod)] # valores das variávies alpha

l_sub = []
l_x   = []
l_I   = []
for i in range(n_prod):
    sub, x, I = gera_subproblem(param, i)
    l_sub.append(sub)
    l_x.append(x)
    l_I.append(I)
    #sub.writeProblem(f"C:\\Users\\x-eco\\OneDrive - ufpr.br\\Documentos\\Alexandre\\DW\\sub_{i}.lp", trans = False, genericnames = False)

i = 0
l_solucoes_I = [[] for i in range(n_prod)] # lista de listas com os pontos usados nos subproblemas: usados no final
l_solucoes_x = [[] for i in range(n_prod)] # lista de listas com os pontos usados nos subproblemas: usados no final 
                  # para, junto dos valores de alpha recuperar os valores de x e I

l_solucoes_x[0].append(initial_sol_x[0])
l_solucoes_x[1].append(initial_sol_x[1])
l_solucoes_I[0].append(initial_sol_I[0])
l_solucoes_I[1].append(initial_sol_I[1])

while True:
    custo_atualizado = 0
    mestre.optimize()
    z = mestre.getObjVal()

    # Duais, fo e variáveis
    pi          = [mestre.getDualsolLinear(c) for c in mestre.getConss()] # coleta as variáveis duais
    var_art_val = mestre.getVal(var_art)

    # Coleta a solução atual das variáveis alpha
    alpha_val   = [[] for i in range(n_prod)] # valores das variávies alpha
    for p in range(n_prod):
        l_val = [mestre.getVal(v) for v in alpha[p]]
        alpha_val[p].append(l_val)

    # Usa os duais obtidos para atualizar as fo dos subproblemas
    l_melhor_sol_x   = []
    l_melhor_sol_I   = []
    melhor_p         = -1
    menor_custo_reduzido = 999
    for p in range(n_prod):
        # Atualiza a fo do subproblema p
        l_sub[p].freeTransform()
        l_sub[p].chgReoptObjective( quicksum(param["Cp"][p][t] * l_x[p][t] for t in range(T)) + quicksum(param["Ci"][p][t] * l_I[p][t] for t in range(T)) -
                               quicksum(pi[t] * l_x[p][t] * param["HMaquina"][p] for t in range(T)) ,
                               sense = "minimize")
        l_sub[p].writeProblem(f"C:\\Users\\x-eco\\OneDrive - ufpr.br\\Documentos\\Alexandre\\DW\\sub_{p}.lp", trans = False, genericnames = False)

        # Resolve o subproblema p
        l_sub[p].optimize()

        # Verifica o custo reduzido (compara com o dual da restrição de convexidade)
        z_sub_p        = l_sub[p].getObjVal()
        custo_reduzido = z_sub_p - pi[T + p]  # considerando a restrição de convexidade na posição T + p
        
        # Se o custo reduzido for < 0, ele trás melhoria
        if custo_reduzido < -0.001:
            # Verifica se ele é o melhor custo reduzido dentre todos os subproblemas
            if custo_reduzido < menor_custo_reduzido:
                menor_custo_reduzido = custo_reduzido
                melhor_p             = p
                l_melhor_sol_x       = [l_sub[p].getVal(l_x[p][t]) for t in range(T)]
                l_melhor_sol_I       = [l_sub[p].getVal(l_I[p][t]) for t in range(T)]

    if menor_custo_reduzido >= -0.001:
        # Nenhum subproblema trouxe melhoria, então a solução é ótima
        break
    else:
        i = i+1
        mestre.freeTransform()
        logger.info("Coluna consegue melhorar a solução do Mestre, adicionando nova coluna...")
        c_nova_col, nova_coluna = gera_col(param, l_melhor_sol_x, l_melhor_sol_I, melhor_p)
        mestre, alpha = add_col_mestre(mestre, alpha, c_nova_col, nova_coluna,melhor_p)
        l_solucoes_x[melhor_p].append(l_melhor_sol_x)
        l_solucoes_I[melhor_p].append(l_melhor_sol_I)
        mestre.writeProblem(f"C:\\Users\\x-eco\\OneDrive - ufpr.br\\Documentos\\Alexandre\\DW\\mestre_{i}.lp", trans = False, genericnames = False)

logger.info(f"Sol. Ótima do Mestre: {alpha_val}")
logger.info(f"Custo: {z}")


# Recupera as soluções dos subproblemas
x_sol, I_sol = recupera_solucao(alpha_val, l_solucoes_x, l_solucoes_I)
for p in range(n_prod):
    logger.info(f"Produto {p}:")
    logger.info(f"  x: {x_sol[p]}")
    logger.info(f"  I: {I_sol[p]}")